{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr1 = np.random.normal(0,0.4,(1000,1))\n",
    "arr2 = np.random.normal(0,0.4,())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = pd.DataFrame(np.hstack([arr1,arr2]),columns=[\"value1\",\"key\"])\n",
    "df2 =\n",
    "df2.columns=[\"key\",\"value2\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.merge(df1,df2,on=\"key\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 读取两个df\n",
    "# 确定两个df的keys\n",
    "# 合并，计算标准分\n",
    "[ \"\".join(i) for i in df1.values.astype(str)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1-0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class VarRelation (object):\n",
    "    def __init__(self,df_x,df_y,x_data_col,y_data_col,x_key,y_key,scale_x=False) -> None:\n",
    "        df_y[\"y_score\"] = self.scale(df_y,y_data_col)\n",
    "        df_x[\"x_level\"] = self.scale(df_x,x_data_col) if scale_x==True else df_x[x_data_col]\n",
    "\n",
    "        df = pd.merge(df_x,df_y,left_on=x_key,right_on=y_key)\n",
    "        self.df = df\n",
    "\n",
    "    def scale(self,df,data_col):\n",
    "        df[\"sum\"] = df[data_col].sum(axis=1)\n",
    "        df[\"sum_standarded\"] = (df[\"sum\"]-df[\"sum\"].mean())/df[\"sum\"].std(ddof=0)\n",
    "        return df[\"sum_standarded\"]\n",
    "\n",
    "    def box_plot(self,xticks,y_ticks,path=None):\n",
    "        plt.figure(figsize=(15,10))\n",
    "        sns.boxplot(data=self.df,x=\"x_level\",y=\"y_score\")\n",
    "        plt.xticks(xticks)\n",
    "        plt.yticks(y_ticks)\n",
    "        if path:\n",
    "            plt.savefig(path,dpi=600)\n",
    "        \n",
    "\n",
    "v = VarRelation(df_x=df2,df_y=df1,x_data_col=[\"value2\"],y_data_col=[\"value1\"],x_key=\"key\",y_key=\"key\")\n",
    "data = v.df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "plt.figure(figsize=(15,10))\n",
    "sns.boxplot(data=data,x=\"x_level\",y=\"y_score\")\n",
    "plt.xticks(\"x\")\n",
    "plt.yticks(\"y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.random.normal(10,1,(1,3))\n",
    "b = np.random.uniform(-1,1,1)\n",
    "x = np.random.normal(0,1,(1000,3))\n",
    "y = np.random.normal(w*x+b,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.68783943,  2.77984315,  1.88882261],\n",
       "       [ 0.12453843,  0.59141656, -2.05416702],\n",
       "       [-0.70202714, -0.67135011,  2.4499405 ],\n",
       "       ...,\n",
       "       [ 0.40843656, -0.36413188,  1.06962749],\n",
       "       [-0.00841013,  0.90588649,  7.60923938],\n",
       "       [-0.57822056,  1.00804956,  1.13859016]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/4999 [00:00<?, ?it/s]C:\\Users\\Public\\Documents\\Wondershare\\CreatorTemp\\ipykernel_7280\\815288189.py:29: RuntimeWarning: overflow encountered in exp\n",
      "  alpha[t,i] = min(1, np.exp((dist_prob(w_[:3],w_[-2],w_[-1])) -\n",
      "100%|██████████| 4999/4999 [00:11<00:00, 437.47it/s]\n"
     ]
    }
   ],
   "source": [
    "def dist_prob(w_,b_,sig):\n",
    "    log_p = stats.norm.logpdf(y,w_*x+b_,scale=sig).sum()\n",
    "    return log_p\n",
    "def replace(alpha,n_value,old_value):\n",
    "    if np.random.binomial(1,alpha)==1:\n",
    "        return n_value\n",
    "    else:\n",
    "        return old_value\n",
    "\n",
    "T = 5000\n",
    "w_b = np.ones((T,5))*0.5\n",
    "sigma = 1\n",
    "t = 0\n",
    "alpha = np.zeros((T,5))\n",
    "prob = np.zeros(T)\n",
    "for t in tqdm(range(1,5000)):\n",
    "    w_ = w_b[t-1,:].copy() # 将斜率、截距和残差放入\n",
    "    \n",
    "\n",
    "    for i in range(5):\n",
    "        if i <4:\n",
    "            w1 = stats.norm.rvs(loc=w_[i], scale=alpha[t-100:t,i].mean()*0.1 if t>100 else 0.5) \n",
    "        if i==4:\n",
    "            \n",
    "            w1 = stats.norm.rvs(loc=w_[i], scale=0.05) # 残差也用正态分布作为先验，但小于0的不计算\n",
    "            if w1 <0:\n",
    "                w1=w_[-1]\n",
    "        w_[i] = w1\n",
    "        alpha[t,i] = min(1, np.exp((dist_prob(w_[:3],w_[-2],w_[-1])) - \n",
    "                        dist_prob(w_b[t-1 ,:3],w_b[t-1,-2],w_b[t-1,-1])))\n",
    "        w_[i] = replace(alpha[t,i],w1,w_b[t-1,i])\n",
    "    w_b[t,:] = w_\n",
    "    # prob[t] =  dist_prob(w_b[t ,:3],w_b[t,-2],w_[-1]))\n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQAUlEQVR4nO3dfYxld13H8ffHbhstVCnu8GDbdYuBJkAgNCMCVQQqyUKBasIfbYRUrNlIAhajYJGE/tsoUTQ+kA2shVBLDBQkPNmGBxsFiru1pVuWh4q1LEV3a408xlr9+sfcleF2Zu6595x77/x23q9kMnPP4/ee+e1nzp57z/emqpAkteeHll2AJGk2BrgkNcoAl6RGGeCS1CgDXJIatWuRO9u9e3ft3bt3kbuUpOYdPnz4/qpaGZ++0ADfu3cvhw4dWuQuJal5Sf5lo+leQpGkRhngktQoA1ySGmWAS1KjDHBJapQBLkmNmhjgSQ4mOZ7kyNj01yb5UpK7kvze/EqUJG2kyxn4dcC+9ROSPB+4FHhaVT0FeMvwpUmStjIxwKvqFuCBscmvBq6tqv8aLXN8DrVJkrYw6zXwJwE/l+TWJH+b5Kc3WzDJ/iSHkhw6ceLEjLvTEPZe/eFllyBpQLMG+C7gbOBZwOuBv0qSjRasqgNVtVpVqysrD7uVX5I0o1kD/BhwY635HPC/wO7hypIkTTJrgH8AeAFAkicBZwD3D1STJKmDid0Ik9wAPA/YneQYcA1wEDg4emvhg8AV5acjS9JCTQzwqrp8k1mvGLgWSdIUvBNTkhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSoiQGe5GCS46NP3xmf99tJKomfhylJC9blDPw6YN/4xCTnAS8E7h24JklSBxMDvKpuAR7YYNYfAm8A/CxMSVqCma6BJ3kZ8PWqumPgeiRJHU0d4EnOBN4EvLnj8vuTHEpy6MSJE9PurrO9V394Ltucx3Zb2b+k7W2WM/CfAs4H7khyD3AucFuSx220cFUdqKrVqlpdWVmZvVJJ0g/YNe0KVXUn8JiTj0chvlpV9w9YlyRpgi5vI7wB+AxwQZJjSa6cf1mSpEkmnoFX1eUT5u8drBpJUmfeiSlJjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAHeOBteSTuXAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqVJePVDuY5HiSI+um/X6SLyb5fJL3J3nUXKuUJD1MlzPw64B9Y9NuBp5aVU8Dvgy8ceC6JEkTTAzwqroFeGBs2k1V9dDo4WeBc+dQmyRpC0NcA/9V4KObzUyyP8mhJIdOnDgxwO52hkX2N7GXitSmXgGe5E3AQ8D1my1TVQeqarWqVldWVvrsTpK0zq5ZV0xyBfAS4OKqquFKkiR1MVOAJ9kH/A7w81X13WFLkiR10eVthDcAnwEuSHIsyZXAnwBnATcnuT3J2+ZcpyRpzMQz8Kq6fIPJ75hDLZKkKXgnpiQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjdpxAb736g8P1rzJJlBtG3IsSMuw4wJckk4VBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUV0+Uu1gkuNJjqyb9ugkNyf5yuj72fMtU5I0rssZ+HXAvrFpVwMfr6onAh8fPZYkLdDEAK+qW4AHxiZfCrxz9PM7gV8ctixJ0iSzXgN/bFV9A2D0/TGbLZhkf5JDSQ6dOHFixt1tbKtGRDutUdE8nuuijuFO+j1JQ5r7i5hVdaCqVqtqdWVlZd67k6QdY9YA/7ckjwcYfT8+XEmSpC5mDfAPAleMfr4C+OthypEkddXlbYQ3AJ8BLkhyLMmVwLXAC5N8BXjh6LEkaYF2TVqgqi7fZNbFA9ciSZqCd2JKUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktSoUyrAp22+1HX5vs2WNtrPvBs4zetYjC+7zKZhLTfBarl2bR+nVIBL0k5igEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIa1SvAk/xmkruSHElyQ5IfHqowSdLWZg7wJOcAvwGsVtVTgdOAy4YqTJK0tb6XUHYBP5JkF3AmcF//kiRJXUz8UOPNVNXXk7wFuBf4HnBTVd00vlyS/cB+gD179sy6u/9v/nPPtZdM3aTppHuuvWTitmeta7P1uzbL2mrf49s4uewsDZG6HI9ZtjXLupP2v1GtXdabtz7jZTtsX6eOPpdQzgYuBc4HfgJ4RJJXjC9XVQeqarWqVldWVmavVJL0A/pcQvkF4J+r6kRV/TdwI/CcYcqSJE3SJ8DvBZ6V5MwkAS4Gjg5TliRpkpkDvKpuBd4L3AbcOdrWgYHqkiRNMPOLmABVdQ1wzUC1SJKm4J2YktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1qtf7wJehTwOlofY7TROmIZbrY6t9bDRvkcd31qZkm82fpTnWVst1aSbVZdnWmlO1Vu9O5hm4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1qleAJ3lUkvcm+WKSo0mePVRhkqSt9b2V/o+Aj1XVy5OcAZw5QE2SpA5mDvAkPwo8F/gVgKp6EHhwmLIkSZP0OQN/AnAC+IskTwcOA1dV1XfWL5RkP7AfYM+ePT12N6zNGjmdbOAzr6ZOQ2x3GQ29pt3nZg2R5ln7EE2Y5l3fZrVNOl5DN5bq0vxL21+fa+C7gAuBP6+qZwDfAa4eX6iqDlTValWtrqys9NidJGm9PgF+DDhWVbeOHr+XtUCXJC3AzAFeVf8KfC3JBaNJFwNfGKQqSdJEfd+F8lrg+tE7UL4KvKp/SZKkLnoFeFXdDqwOU4okaRreiSlJjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqP63shzyunTzGjSul22vYxGVZNMW/csTZJmXX/WYzo+bVJzs1n33cV2/J2fik7FBl6egUtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqN6B3iS05L8Y5IPDVGQJKmbIc7ArwKODrAdSdIUegV4knOBS4C3D1OOJKmrvs2s3gq8AThrswWS7Af2A+zZs6fn7rYHmw/NV9+mYEM2LTq5ry7b67LsEA3PNlt+Uo2bHZcu25j2mPZtbjbNfro872WaZgxNa+Yz8CQvAY5X1eGtlquqA1W1WlWrKysrs+5OkjSmzyWUi4CXJbkHeA/wgiTvHqQqSdJEMwd4Vb2xqs6tqr3AZcAnquoVg1UmSdqS7wOXpEYN8ok8VfUp4FNDbEuS1I1n4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNWqQ94Ev00aNajZrXrPspjbztl2e33gdy6hrliZUQ+xvmm3O2rhqq4ZTW82fZh/TzDu5v1mf8/p6Z2lENmSzqFmads2zWdcknoFLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1Kj+nwq/XlJPpnkaJK7klw1ZGGSpK31uZX+IeC3quq2JGcBh5PcXFVfGKg2SdIW+nwq/Teq6rbRz98CjgLnDFWYJGlrgzSzSrIXeAZw6wbz9gP7Afbs2TPE7k5J26UR1dCW/byW2UhrXtucpmnUPddesum6k7Yz9PPo2mRu/eNJjbsmTdtq/Xk2BdusnqH1fhEzySOB9wGvq6pvjs+vqgNVtVpVqysrK313J0ka6RXgSU5nLbyvr6obhylJktRFn3ehBHgHcLSq/mC4kiRJXfQ5A78IeCXwgiS3j75ePFBdkqQJZn4Rs6r+DsiAtUiSpuCdmJLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNWqQZlaShjdEM6RFNFSaZ1OsjZpyTdNkaqtGWdPUN77M+ho2ar61qCZqnoFLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1Kj+n6o8b4kX0pyd5KrhypKkjRZnw81Pg34U+BFwJOBy5M8eajCJElb63MG/kzg7qr6alU9CLwHuHSYsiRJk6SqZlsxeTmwr6p+bfT4lcDPVNVrxpbbD+wfPbwA+NLs5c7VbuD+ZRexBevrb7vXaH39nMr1/WRVrYxP7NONcKNPpH/YX4OqOgAc6LGfhUhyqKpWl13HZqyvv+1eo/X1sxPr63MJ5Rhw3rrH5wL39StHktRVnwD/B+CJSc5PcgZwGfDBYcqSJE0y8yWUqnooyWuAvwFOAw5W1V2DVbZ42/0yj/X1t91rtL5+dlx9M7+IKUlaLu/ElKRGGeCS1KgdEeCTbvlP8stJPj/6+nSSp6+bd0+SO5PcnuTQkup7XpL/HNVwe5I3d113QfW9fl1tR5L8T5JHj+Yt4vgdTHI8yZFN5ifJH4/q/3ySC7s+twXVt+zxN6m+ZY+/SfUte/ydl+STSY4muSvJVRssM58xWFWn9BdrL7D+E/AE4AzgDuDJY8s8Bzh79POLgFvXzbsH2L3k+p4HfGiWdRdR39jyLwU+sajjN9rHc4ELgSObzH8x8FHW7l141snf7yKOX8f6ljb+Ota3tPHXpb5tMP4eD1w4+vks4Msb/BueyxjcCWfgE2/5r6pPV9V/jB5+lrX3tG+b+ua07rzquxy4YeAatlRVtwAPbLHIpcC7as1ngUcleTwLagcxqb4lj78ux28z2+L4jVnG+PtGVd02+vlbwFHgnLHF5jIGd0KAnwN8bd3jYzz84K53JWt/KU8q4KYkh0dtAZZV37OT3JHko0meMuW6i6iPJGcC+4D3rZs87+PXxWbPYRHHb1qLHn9dLWv8dbYdxl+SvcAzgFvHZs1lDPa5lb4VnW75B0jyfNb+Af3suskXVdV9SR4D3Jzki6MzgkXWdxtrvRC+neTFwAeAJ3Zct69p9vFS4O+rav3Z0ryPXxebPYdFHL/OljT+uljm+JvGUsdfkkey9sfjdVX1zfHZG6zSewzuhDPwTrf8J3ka8Hbg0qr695PTq+q+0ffjwPtZ+y/PQuurqm9W1bdHP38EOD3J7i7rLqK+dS5j7L+vCzh+XWz2HLZNO4gljr+Jljz+prG08ZfkdNbC+/qqunGDReYzBud5cX87fLH2v4yvAufz/RcJnjK2zB7gbuA5Y9MfAZy17udPs9aBcdH1PY7v33T1TOBe1v5yT1x3EfWNlvsx1q5TPmKRx2/dvvay+Ytwl/CDLyB9bprntoD6ljb+Ota3tPHXpb5lj7/RsXgX8NYtlpnLGDzlL6HUJrf8J/n10fy3AW8Gfhz4syQAD9Va17DHAu8fTdsF/GVVfWwJ9b0ceHWSh4DvAZfV2m9/7u0MOtYH8EvATVX1nXWrz/34ASS5gbV3SuxOcgy4Bjh9XX0fYe1dAHcD3wVetdVzW0J9Sxt/Hetb2vjrWB8scfwBFwGvBO5Mcvto2u+y9od5rmPQW+klqVE74Rq4JJ2SDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUqP8DcIlJJIlwoegAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(stats.truncnorm.rvs(a=0.1,b=2,loc=0,scale=1,size=1000),bins=200)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3033287274197666"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha[:,4].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqS0lEQVR4nO3dd3wUZf4H8M83nRJCSygJEAKhSRNCLyII0jzshwUbWH6e7fQKnoKn6Ilyet7hKWI5PTzFLiqISpcmBAgQOoQIoSUBQguEJPv8/tjZ2dnd2ewm2bCZzef9evHK7szs7DPL7Gee55lnZkUpBSIisr6wYBeAiIgCg4FORBQiGOhERCGCgU5EFCIY6EREISIiWG/cuHFjlZycHKy3JyKypA0bNuQrpeLN5gUt0JOTk5Genh6stycisiQR+dXbPHa5EBGFCAY6EVGIYKATEYUIBjoRUYhgoBMRhQgGOhFRiGCgExGFCMsH+tqs49ibezbYxSAiCrqgXVgUKONnrwUAZE8fE+SSEBEFl+Vr6EREZMdAJyIKEQx0IqIQwUAnIgoRDHQiohDBQCciChEMdCKiEMFAJyIKEQx0IqIQwUAnIgoRDHQiohDBQCciChEMdCKiEMFAJyIKEQx0IqIQwUAnIgoRDHQiohDBQCciChEMdCKiEMFAJyIKEQx0IqIQwUAnIgoRlg50pVSwi0BEVG1YPNCDXQIiourD2oEe7AIQEVUj1g50QxWd3S9EVNNZOtBtyvwxEVFNZOlAV4ZOFxtr6ERUw/kMdBF5T0RyRSTTy/zbRGSL9m+1iHQLfDHNKZcaOgOdiGo2f2ro7wMYWcb8/QCuUEp1BTANwOwAlMsvxgxnnhNRTRfhawGl1AoRSS5j/mrD07UAkgJQLr+cKSrWH7OGTkQ1XaD70CcC+D7A6/SquNTYh36p3pWIqHryWUP3l4hcCXugDyxjmfsA3AcALVu2rPR7GocqsoZORDVdQGroItIVwDsAximljntbTik1WymVppRKi4+Pr/T7GjP8XFFJpddHRGRllQ50EWkJ4EsAE5RSuytfJP8Za+UXim2X8q2JiKodn10uIvIxgCEAGotIDoBnAEQCgFJqFoCpABoBeENEAKBEKZVWVQU2MvabF5cy0ImoZvNnlMstPuZPAjApYCUqB2MNnYFORDWdta8UdQl0nhQloprN0oFu7HI5XHA+eAUhIqoGLB7ozkTnKBciqumsHeiGbnN2uBBRTWftQOf90ImIdJYOdN6ci4jIybKBXlJqwz8X79Gf814uRFTTWTbQ52UcxqIdx/TnvJcLEdV0lg30ohLXC4nYh05ENZ1lAz1MXJ8zzomoprNsoItboJeyE52IajjrBjpcE509LkRU01k20MEuFyIiF5YN9DBxr6Ez0omoZrNsoIvvRYiIahTrBrp7lwsr6ERUw4VOoLMXnYhqOMsGumcfepAKQkRUTVg20N0xz4moprNsoAtr6ERELiwb6O7DFNmHTkQ1nWUD/dWfdrs8Zw2diGo6ywb6r8cLg10EIqJqxbKB7o5XihJRTRdCgR7sEhARBZdlA71x3SiX58xzIqrpLBvoDw9NdXnOGjoR1XSWDXTPXyxiohNRzWbZQHePb9bQiaims26guwU48zx4ft6Th3X7TwS7GJaklMKF4tJgF4NChIUD3T3RGenBMuHddbj5rTXBLoYlzfhhF7o9+yMultiCXZSA2HX0DH7JOh7sYtRYlg30/LMX9ccigauhn75QjF1HzwRobURle2PZPhSV2HD6QnGwixIQV7+2Ar+dvTbYxQio7PxzOHHuosf03cfO4PjZoiCUyDvLBvrrS/fqj8NEAlZBv2X2Wlz92opyv25t1nEkT56PbYdPlfu1pTaF4lLr1NB+3pOHc0UlHtPLCqVfso5jX97ZqiyWpQW6gXn01AUkT56PrzblBHbF1dSi7ceQPHk+DhecD/i6h/x9GXpM+8lj+oh/rMDwf5Q/K6qSZQPdSADYAvCNKLUpbDt82mVacakNR09d8PnamUv2AADW7HM2N1fvzcc976+HzaaQnn0CJV5C+5731yP1qe8rUXLfvtqUg283H8Z/12Rj1vJ9FV7P3tyzmPDuOvzD7V46APDHzzajqMS1P7ik1IY9x87gt7PXYtgry13mldoU9uYGtzV0rqgECzOPBLUMgO8rnX/cdhRvlfH/9uO2o8g97dxPs7SD5yfrD5ouv3Rnrl/nPbbkFCB58nws3nHMdP7mgwV4Zl5m0K/UnqttZ+ah8leoKsOs5h5MoRHoAtgUcORU+Y/Ow19djt4vLMKOI6fR5i8LXOY9++02pD71Pfq+uBgFhd7/42w2hVV77UFu3K/vn7MBS3bmYsWePNw4aw36vrgYAJB3pghz1v6KJz7djLVZx7F8d165y11ev/9kMx7+eBOmztuG6d/vrPB6HJ/DxgMnPeb9sO0YbnhzNX7a7vzyz/hhl9dazNiZK3HVqyuQc9LzvjxfbMhB9+d+dGnSrtyT7zVYACB58nw8+eVWv7cFAKbMy8QDH270GQQHjhfig9XZ+HDtr6bzS23+BVrhRc+WDWDffwHgVKF5K+e+ORvwotv/W1FJKYa9sgzLd+fhvjkbMH72WhSVlOLjdQewWqtYFJeal+vu99fj5rfW4KQhkI6dvoC3V2ThjKGl9ZvXVwEAXlu0x3Q9f/p8Cz5Y8ys2HijAc99u16enZ5/A+Nlr9ANLWfbnn8NDH2003Q+Mlu7MxfGzRZj4/nqvFYH75mzw+X7lsd1QwTtUydr/qfPFfn0elREagQ7BWyv2od+LS3x+YAdPFOq1yC05BdiTexa5Z4rw6/FzHsv+Z1W2/vifi8136Fd/3OXS/bMu+wS+33oEu4+d0X/J+ohWw88/exE9p/2EK2YsxZSvM/HFxhw8/kmGX9t406zV+M+q/X4tW2pTsPkImOz8c0iePB/r9p9A3pkifaSFzaaw7fAprzUuR3CFu18IoMk8dBr3/jcdX286hLnrDuCtFVley7DjiP3LUmIIndMXinH/nHQ88dlmFBQW49hpZ6Df/u4vmPhBul42pRSSJ8/Hc99u10Po43UHvL7f6n35eNtQnrNFJfhy4yEA9i+bNzabwuAZS/HMN9vw9NeZ+vSMgwX4eN0BHCo4j9SnFuCdn7Pw76V7MeHdX0zXM/zV5eg09QfszbXvo8banU0pfLv5MLo99yNW7skHYK9tTvpgvctB7ec9zoN/zsnz2Jd3Dn/RDmJZ+efw6fqDePLLrfo+ueHXkzh/sRQvLdyJJz71bEE9qu1/Sim8tmg3XliwA0t3eVYwakWG649zT1/AoJeXYM2+43oXxw1vrsZ7hv3zxllrsDbrhMt3yJuP1x3Ad1uOYNF284O1zabw5cYc3P3+evR8fhEW78zFlK+3YW/uWSzMPArA8ycpA+Wa11fqjwdMX+Lx+RktzDyKOV4O+ADw2NxNGPrK8io9AR5RZWu+lMRZM845eR4p8XWx/fBpHDxZiKsva6ovVlB4EYNeXoobeybhjn6t8OD/Nurz3LtV3GsL/1mVjSljOuF8cSnqREcg89Ap3D9ng8dR+6ftx1xqqID9S+Vw3K2JdtjwvheKSxETGQ6lFPLOFuGbjMN4fv4O7Hp+JNZnn8T67JOY0LcV/vfLARwuOI8nR3cEYK+5vLZ4Dz6a1Ad1oiPQ64VFiIkIw+onh3n9yFbts4fGV5sO4ZP1BzC6SzO8fmsPLN6Zi3v/m45HhqUiOiIMv+3VAo3rRgMAlu/Ow5w12QCA9dmeNXSjx7wcqHYfO4N2TWJdDjgK9q6tc0UleHHBTvywzfn5fbUpB2eLmqJ364b6tJcW7sLkUR30g8t7q/Zj2a5cfX7moVPonBgHpRRW7MnHoLaNkZFTgFvftgftxIGtERYmmLPG+eVznMNYmHkU9WpFoH+bxigovIgzF0qQUC/aYzveXpGFFxbsAAB8+WB/2BTw/PwdHsudKyrBvxbvwe+Ht8MeLchzThbiqa+24trLE/XlbErp+8nuY2cQGxOBcf+2145/95H5fiqG1zp8vM6zi+X91dl4c5m9uyYtuQEiDAfjFbvz8PinGSgutR9QAPPun3XZzu6ZTQcLcPDEeXsFw0eQ+nNuKONAAQAgLEzw9Ndb0b5pPQztkIDoiDA0rhuNd1Zm4W8LXFsnCgqTPliP7OOF2Dx1hMd3zh+5py/g/g834K3beyKhXoxpud1bXtn5hUiJr+MSyjabggjwwIf21sGEvq2c5VQKRSU2xESG6wfKdk9/j+zpY8pdXn+ERKAb9ynHxz/6Xz8DAHY/PwqXPbMQA9s2xnPjOgMAPt+Qg883uJ4semOZa/9k5iHXvnQAeGHBDry7cj92ThuJsTNXesz3xv29vJm/5Qhu6JmEl3/YpX8BAeDMBWcz/b1V+/Wd+8nRHXH6QjHufn89AOCyZ37AoscH6zW/JTuP4cr2CR6/7gQAF4oNO6QCvttyBK/fCpzWaqr/0lokjnB8cEgb3PneOr+32ZvfvL4SO6eNwiJD10mpzYbH5mZg/tYjGNu1mcvyb/+8H2//vN/lCzBr+T5MHtUBJYYvW1a+s4U1duZKZEwdjl/2n8D9czagcd1o5Btqud9tPYIxXZqZ1uocX8rs6WNwxYxlOHW+GB9N6uOxnCPMAWD2cs9WiFIKIoI3lu3FWyuyUK9WpD7v8w05+GX/Cfxi6MNWyhmkIvZuCIe1Wc7lVu87jqW7cvHGbT1RpIWKsVvF7FzSSwudYWjWJeVopRjXsXpfPj5Yne2x7IXiUr1FlX+2yFeeQ2ldoev2n8CjczMwaWBrPD22kz4/89Ap/WAxdd42ffoUwzquMxz4HIyfSbfnfnSZt/Poabzz8350TYrDZc3j0LNVA31eUUkpHpizAQNT41FQeBGbDhSg998WY/+Lo7HxQAHyzxahTlQEbvfSyrr6tRWY0LcVsg0t+hS3rtprZq7E1kOnkD19DN5cvg8vL9yFqYZtrkqhEeiGvcqmFOZlOHfQmUv2oLhUYemuvDL7OXPPuA4/yjvjeSL03ZX2JuUd71Y+2Mw88dlmzF6RhV3HXPsH055fpD/ONXRB/Hr8HK6Yscxl2bv+s15/fM/76fjo3j7o27qRx3tN+87e3+neRREd6doL933mUXyfeRTje7Uo38Z4caHYhnkZh/CWIQRLbArzt9pPTPp7am1exqEyT4R3f+4nTLvWfgDPdxta9sjHm3DwRCEiw507jvuqikpK9W6YW98x/3I7LNx21GNacalCQWERjmvDa2f8sEuf990Wz5Ow7tvirQvhq032ffvdlfv1/8NSm/PgvDMAQ25nLt7rcoB0UEphwPQletk2HihATGTZvbYKCnf/Z71erndW7sc72vdo57SReGTuJp/lcWyzv0a+Zq/MOSpS47o3R3r2SayaPBRbc05h6a48LN2Vh0eGttVfs+3wadzw5moAwN0Dkstcf1ndKgCwVTsf89LCnXrF7LnvtrssY7MphHnptqwMn4EuIu8BGAsgVynV2WS+APgngNEACgHcpZTa6L5cVTLm9N2GQAOAmUuc/dtD/r7M73VOMdQW3Bmbn4HmHubuHF8GAB5hDti7nFzWd/SMy8mqsqQ9v8hr37lZX/jqffno36axX+s2enRuhsu9eLLynOGxZEeuySs8TyY+OjfD5/uEl9Gxmp1/Dp8ZWk6//zQD65+6Sn++3KQfuTxKbDb0n77EpRVRFqWcBzOB7+2bZgiIk15OpJalfZNYr/uaWZgD9hOk7l2GxpaemU/TvbdOc04WuvzfV5V5GfaupLVZx/HnL7bo088WOfvDjS3uQA3YeXOZ91FJq/cdx8DU8n93fPHnpOj7AEaWMX8UgFTt330A3qx8scqJF4l6FR4mftfa8s8WeXxhHd41HEgcbn37F5eTdOURFeHc9f5lOOF83stl8G8sLf9Qy7985X3Ey2du3WAFha4XlJU1WqL1k/N9vnePaT/5HeaAvYbuqKX/1c8DcGW4t8T8sTXAQwJLbUCb+DoBXWdZxs9ei6QGtfTn73kZZFBi893v3y0prlJlSf+1aiqFPv9XlVIrAJT17uMA/FfZrQVQX0SalbF8wNWNCYmeoyrx3zVlNw/95a27akIFu58uFNuQ0tj+ZfZnOFhFLtgqr6fKOAAY+VOD81VzdVeO7A+ImAjnqJXUhLqX9s01JTYbmsbFXNJQT6pf2+cy+3J9txo251Ruf6yqCwkDMWwxEYDx1HqONs2DiNwnIukikp6XF7ix1zf2TArYuqykVSPfO6djiFx15GjaG0/6emM2lC7QHF9SY+uhPKLCK/51Ukrhw7Xeh1wGmmMbB7RthE/u73fJ3tfIZgO2HDyFutGXrkL2Sbr5hVZGay7BvWi8XR9QWYEIdLOOStPSKqVmK6XSlFJp8fHxAXhr4Poeifj9Ve0Csi6reXJUx2AXAQAwqIJ9ga0a1a7wa/3x8o1dK/S6yAqerGpbiZru9iOeo6qqkuP0QmR4mM8Tm75c3rJ+hV5XcP4izhSVXPLWSXVQVWPRAxHoOQCMQyCSABwOwHr98tINXVErKtz3ggC6tagPAHhwSBt9WrCam5X13cMDMbJzU6x5cmiF11Hb5HO7q38y6teONFnau+gI/z5/dymN6+A/d/Wq0Gsdhndq4nXezWkVG5kTaVJD75XcwGRJV83iYnwu440/J3kDydFtFBkehtpREfj43r4VWk9EmODzB/rjmm7NTfenP4xoh53TzE/BObrrbundskLvXVHG0S3+evY3lwW0DMbrJgIpEIH+DYA7xK4vgFNKqUt2cwxHXWrdX1wvormtT0uPGtqV7e2tgr4pzmF8kRVoJteNjsDnD/hupnq7mrKyvvi//uicaD8pUzuqYs3V5X8cgg1PD/eY/siwVCz7w5AyX1vP7ZyFPzW8sV2b4ekx9hbF9dq44oTYGESEh+HOfq0wuktTl+XfuSMNDXwcWO7qn4y370jz+d5mHGUxExHmuT3je7XEE8M9W4KdmtXTt//hYal447YeFSpPRV3bvbnHtM3PjMAV7cpuAZ/TRg05hm72a9MIr97crVzvXScqHCv/PBThYYKZt1yO7c95Bnf3Fg0QE+n9gD9lbCeM7Vb+U259Uxr6XshE1t9GIyW+7Eqce4Vm8RNX4Poepr3IPt3UMwkPXNHGY3pFu/V88blWEfkYwBoA7UUkR0QmisgDIvKAtsgCAFkA9gJ4G8CDVVJS7+UDAMTHRuO+wSn69Beu6+Ix/vqmtBbIfPZqDDbs7Hf2b4XyuvqypoiP9bx60N0tvVvgngGtcXtfew3kmm7N3ea3xENXlq+2kD19jMuFEu4HjWeu6YQdz410ubISAD6c2Acdmsbio3v7YNXkoWjVqI7HWOcuiXFoWCcK9WtHYd7vBuAPI9q51NwSYqPx9JiOmDjQ/jkPaR+Pjs3qoUui6xn/UZ2boney5xdu0qAU7H1hFF65uRtevL4Lpl5jv9ji2XGd8e9bXYOwV3JDbJo6AjunjcT6p67C67dejmEdEvT5G6cMxxQ/LtZo0dA5quHqy5qgdeM6yJg6HJMGpWDyqA4AgI7N6rm8ZninBH1Y5W19WiL96atwg5fzNI3qRukHgHARjO7SDH8Y4Rn8ZgH795vKF6AO3z08UH+c2iTWY35crUh8cE9vZE8fg4e91EZbNbSffzFWaEZ38QzWmMgwxEZHYPETV2DR44P1lkpsTAQ2Th2Opl5aJY7vWLHJiBFjON4zIBmxZfShd2jq3L55vxugP05NcE5vXo6WUViYINEw0sVMuIh+8H7phi5oE1+3QgHcpF40nh7byeN7Nii1Mb58cID5iyrJZ/VOKXWLj/kKwO8CVqJycnzxRAR/Gd0Rsw3jpR2XbIeHCUptCgmx0R418tFdmuHPX5Q9uuG7hwfiyS+3ugzbio+NRq3IcNycloROzevp63jtt931y94fujJV3+EfHpqKRnWi9MurAeDmtCSs2puvP68TFY5zF82H7b1yUzePkAY8x1rHx0ajVlQ4Pr2/Hz5c+6t+75EuSXFY+Nhg19e6HQyMtZ5uLeqjW4v6LiNQXry+C4Z1bIJThcWIigjDnf1boXZUhEvzsVdyA0y/oSviakVi1d583OZ2UU6E9vm7N7NFBFl/G40fth1F9vFCxGm1pJjIcMREhmNs1+YY3qkJ2j+9EADQsE6U6efk7qsHB2D74dP4afsxPDm6g0uL5t5BKRh5WVMkNqil3+1yTNdm+Mvojnjx+q74Jes4uibV17v0HF/qB4e0QduEunj8081IrF8LW7STqY7/ioeGpuLvP9rvRvnNQwNw29u/4JFhbV1uwnZL7xa4sWcSzl4oxl+/3Y4WDWvh4Inz6NSsHtok1MW+3LPYfuQ0ZtzYFX/8fAu+/t0APDZ3E7KPFyI6Igz/+G03rNxzHNdenoidR88grVUDPPON57UT9w5OQduEutiSc8pl6Gmv1g3xdcZhDGjrPIfhXpOe/8hApCbEIjJc9IrTKzd1x+AZS3Fl+wTTrrYh7eOxbFcepo27DC8u2Il+Wmu4W1IcNuecwqzbe6BudCS+3HgI/VIaeVzF/Py1ndGxWT3syz2L2tHh6JpYH4NnLEVK4zp6lynguw/68pb1sUm7pcDu50fhnZVZ+p1QeyU3RMbU4Zj8xVbTi8LCwgQPD0vFvYNT9M/E7IR3Qmw0cs8UoW9KQ5crVx0uax6HuFqercw5Ez2vPA4Uy4/3M7us3SEmMtznPRPCRLDsD0MQVysSh0+dR4PaUeg/fYnLMp0T4/DtwwMxa/k+TP9+J0TsXR07tL5B4xc1zdDXWifaucM3Mdwr4pGhbfH4iPYAnGfUJw5sjSljOyF5susYZ0d/32+6NzftHjL2DtSNjkC3pPr689v7ttIDPcKk+yfM7bNrVNez1ZFYv5Z+6bxjXHVc7Uj8n+E8hHE9M2/poe/EA9o2xnPjLnO5pLssYWGCUSa1RAdfffWxMREeI2Ya143G4HbxLq0yh/AwQbI2dLJ+7UgUFBbj2u6JiI2xl79PimsLb1z3RPx6ohC39W2FxPq1EBMZjivbJ+j3ETF+DpHhguJSha5J9bH12atd1vP+3b30IL2jXzI6NKuHA8cL8acvtqBDs1i8enN3LMw8igc+3IARnZripun2cwGOA0qpUrju8iRcd7m91TDzlstxobjUNNDrxURiXPdEjOueqAd6bHQEbk5rgcGp8S7jsgH7/7fjIH5Zc8+x1i0b1cb/JvVBj5bm5xTemtAThUWlaFAnCrMm9NSnX9OtOTbnnEKXpPrYo13QFBHuuU/ert0HxdgK3frXEfq+v/WvI/DHz7bg8RHt9BErkwaleFyJ+dINXTFCu8tnVEQYHhzSFg8OcbZW6teO8nj/n/90JQa9vFRvXRoPcCKCtyb0RNuEuhg/ey3yzhTpB/CXbuiKa/+9ChHhYcgzXHEubn8BeLRmA83yge7usatSXY7kvhi/1A20Wt9VHRNwZYcEPPVVpsuyDWvb57vvho7ng1IbI6mBcyihWV+s+wHGEbSOdfz1mk74JD0HM27silPni11qUGaiwsPw27QW6NGqPm7q2cLjcuLRXZpiwdaj5gcDw6L/HN8dY7yE6Zu398BNs9a4HCyMHDt27+SGHk3w+rX9q0n7a82TQ3HWyzDH5X+8Egu2HsFri3broeyvJ0a0x5SvMzGwjM+7aVwM/nZdF/25o4virPZjH8bL9zdOGQ6z61OGd2qCIe2dXUdhYYK+KY1w4IT9ZnCbDxYAAEZ2buqxrzgOGGbrNTtgu5t2bWfsyz2LKWM7ITxM0KKh57BXf+5aWNY+GR0RbnrgnTiwNW5Ka4G4WpHYod2S1ljmLolx6GXSTQfA5f8yNiZSP1Dc2a8VBreLx7COTfBp+kH9AjrH5zbjxq6mNWcHYwv1gSvaoEXD2pj/yEC08dLH7rjR35yJvXHo5HlM1u6LEx0Rjk1TR6CopBTZ+YU4d7EE17+x2mPUU0rjOvjMj3NvlRGCgV6+IYxmO/A7d9pHXrgHeqO69nByHAAcHEd69x3ZJM89OGrujr93DWiNuwa09v1CjYjgpTKG5714XVdMGpRi2gdobN2M6+79pE+v5IZltnTaN43FFe3i8ehVqR7z2mpfjisNIVYZzeJqAV4qOQ3rROH2vq30Wl55TOjbCrf3aVlmi8+b+wanYOaSvYg2fMZmB5T9L472un5Hjbes8cl1tL5mZTIq2J8T8BP8+FziY6M9bh8RCCKit9wcLb0IQyXjW8N5AX89O87jTiQubkprgZvKGOn01JiOiAgLwwvXddZr42atEncdmtZDh6b1oJQ90B0ffXREONprff7/m9RH7yJ13K104qDWZZ4gDoSQC/Tycu92MHrzth4ufchD2ifgn+O7ewyV69u6ER4Zloo7+tm/MF0S47D10Cm/LjT5TbfmiI4Iw/BOTX0uWxFxtSO9No8BoEfL+mWGuT8SYmPwwT29Ted1al4PG56+yu8+74q4Z0BrpATgasOKhDkAPDS0LYa0TzA9Qenv+h2fT1kXi71+6+WYu+4gOrmdxPW17vJ4a0JP9H5hcUDW5Y2jJRNp0uVSUY7Pb0Y5rj1IiI3BK+Uc2WN0a5+W+NfiPahtclLX2Iq5s38y4mpFmt41MtAkWD8dlZaWptLT0yv8ekdfc0XvK+x4/d4XRrnUFALh1PliHD11QT9aE/nj5z156JpYXz8hXF7Jk+ejQe1IbJo6olLlqOx3y5eiklI8M28bHh/ezuM+5BWVf7YI3289ggn9kgOyPn/YbAonCi/qNfBLRUQ2KKVMx+vW2Br6R5P6YF7G4SoZKx5XK9L07DZRWQalVu7q6bfvSEOn5p6194owO4kcKNER4Zh+Q8Wu4vWmcd3oSxrmgP38x6UOc19qbKD3b9sY/X2ccCSykrKumi2PzVNH+H31NVUvNTbQichcRbt8KPgs+yPRDWpHVtnls0REVmTZGnqrRnUQy/ugExHpLFvFVQjcUC0iolBg2UCHUj5/cZyIqCaxbKDba+jBLgURUfVh3UBX5j+VRERUU1k20LceOoWC88XBLgYRUbVhyUAv1H5txXG/YyIismigV9UvZhMRWZklA93k7qFERDWeJQPdcT/osn6LkIioprFmoGs19MdNfoyXiKimsmaga385bJGIyMmaga5V0XnpPxGRkzUDXfvLPCcicrJmoGuJzjwnInKyZKDrWEUnItJZMtAVB6ITEXmwZKCDXS5ERB4sGeg8KUpE5Mmaga7X0JnoREQO1gx0OMahB7kgRETViDUDnX3oREQerBno2l/W0ImInKwZ6I5L/1lHJyLSWTTQtQfMcyIinSUD3YF5TkTkZMlA10+KshOdiEhnzUB3DFsMcjmIiKoTSwa6AyvoREROfgW6iIwUkV0isldEJpvMjxORb0Vks4hsE5G7A19UJ8V7cxERefAZ6CISDuDfAEYB6ATgFhHp5LbY7wBsV0p1AzAEwCsiEhXgsuo4Dp2IyJM/NfTeAPYqpbKUUhcBzAUwzm0ZBSBW7Gcp6wI4AaAkoCU1vhnHoRMRefAn0BMBHDQ8z9GmGb0OoCOAwwC2AnhUKWVzX5GI3Cci6SKSnpeXV8Eis4ZORGTGn0A3i033XuyrAWQAaA6gO4DXRaSex4uUmq2USlNKpcXHx5ezqMb1VPilREQhy59AzwHQwvA8CfaauNHdAL5UdnsB7AfQITBFNOO42yKr6EREDv4E+noAqSLSWjvROR7AN27LHAAwDABEpAmA9gCyAllQI95tkYjIU4SvBZRSJSLyEIAfAIQDeE8ptU1EHtDmzwIwDcD7IrIV9pz9s1Iqv6oKzT50IiJPPgMdAJRSCwAscJs2y/D4MIARgS1aWeWx/+UoFyIiJ0teKcpfLCIi8mTNQGcfOhGRB2sHOhOdiEhnzUAHf+GCiMidJQPdgTV0IiInSwY6rxQlIvJkyUB3YAWdiMjJkoHOn6AjIvJkzUDnT9AREXmwZqBz2CIRkQdrBrr2l4FORORkzUDnLxYREXmwZqA7HjDPiYh01gx03suFiMiDJQOdv1hEROTJkoHOGjoRkSdLBvry3XkAOMqFiMjIkoE+c8leABzlQkRkZMlAdwgPY6ATETlYOtAjwxnoREQOlg70Ehvvo0tE5GDpQLcx0ImIdJYO9FL+0gURkc6SgR5XKxIA0DelUZBLQkRUfUQEuwAVkRAbjQFtGyEy3JLHIyKiKmHJRCy1KYSHWbLoRERVxpKpWGJTiOAYdCIiF5YM9FKbQhiv+ycicmHJQC+8WMIaOhGRG8sF+rmiEpwsLMbFUluwi0JEVK1YMtABoGXD2kEuCRFR9WK5QHdcHNo0Lia4BSEiqmYsGOiOH4gmIiIjywW642J/jnIhInJluUB33JCLeU5E5Mpyge64Hxdr6EREriwX6HofOvOciMiF5QKdfehEROb8CnQRGSkiu0Rkr4hM9rLMEBHJEJFtIrI8sMV0Yg2diMicz9vnikg4gH8DGA4gB8B6EflGKbXdsEx9AG8AGKmUOiAiCVVUXigt0FlDJyJy5U8NvTeAvUqpLKXURQBzAYxzW+ZWAF8qpQ4AgFIqN7DFdHJcWMQ8JyJy5U+gJwI4aHieo00zageggYgsE5ENInKH2YpE5D4RSReR9Ly8vAoVmKNciIjM+RPoZsnp/mOeEQB6AhgD4GoAU0SknceLlJqtlEpTSqXFx8eXu7CAsw+dN1skInLlz0/Q5QBoYXieBOCwyTL5SqlzAM6JyAoA3QDsDkgpDZwnRZnoRERG/tTQ1wNIFZHWIhIFYDyAb9yWmQdgkIhEiEhtAH0A7AhsUe0cXS6McyIiVz5r6EqpEhF5CMAPAMIBvKeU2iYiD2jzZymldojIQgBbANgAvKOUyqyKArMPnYjInD9dLlBKLQCwwG3aLLfnMwDMCFzRzOl96Ja7JIqIqGpZLhbZh05EZM6CgW7/yzgnInJluUB3jJhkHzoRkSvLBbqNJ0WJiExZL9BtvLCIiMiM9QLdcY0qA52IyIXlAl2xD52IyJT1Ap196EREpiwX6Lw5FxGROQsGuv0vK+hERK4sF+iKV4oSEZmyYKDb/7IPnYjIleUCnX3oRETmLBjo9r/CgehERC4sF+jOPvQgF4SIqJqxXKDzXi5EROYsF+iKP3BBRGTKcrHIPnQiInOWC/SmcTEY06UZYmP8+vU8IqIaw3Kp2LNVA/Rs1SDYxSAiqnYsV0MnIiJzDHQiohDBQCciChEMdCKiEMFAJyIKEQx0IqIQwUAnIgoRDHQiohAhjnujXPI3FskD8GsFX94YQH4Ai2MF3OaagdtcM1Rmm1sppeLNZgQt0CtDRNKVUmnBLselxG2uGbjNNUNVbTO7XIiIQgQDnYgoRFg10GcHuwBBwG2uGbjNNUOVbLMl+9CJiMiTVWvoRETkhoFORBQiLBfoIjJSRHaJyF4RmRzs8lSGiLwnIrkikmmY1lBEfhKRPdrfBoZ5T2rbvUtErjZM7ykiW7V5/xKpnr+gLSItRGSpiOwQkW0i8qg2PZS3OUZE1onIZm2bn9Wmh+w2O4hIuIhsEpHvtOchvc0ikq2VNUNE0rVpl3ablVKW+QcgHMA+ACkAogBsBtAp2OWqxPYMBtADQKZh2ssAJmuPJwN4SXvcSdveaACttc8hXJu3DkA/AALgewCjgr1tXra3GYAe2uNYALu17QrlbRYAdbXHkQB+AdA3lLfZsO2PA/gIwHehvm9rZc0G0Nht2iXdZqvV0HsD2KuUylJKXQQwF8C4IJepwpRSKwCccJs8DsAH2uMPAFxrmD5XKVWklNoPYC+A3iLSDEA9pdQaZd8b/mt4TbWilDqilNqoPT4DYAeARIT2Niul1FntaaT2TyGEtxkARCQJwBgA7xgmh/Q2e3FJt9lqgZ4I4KDheY42LZQ0UUodAewBCCBBm+5t2xO1x+7TqzURSQZwOew11pDeZq3rIQNALoCflFIhv80AXgPwJwA2w7RQ32YF4EcR2SAi92nTLuk2W+1Hos36kmrKuEtv2265z0RE6gL4AsBjSqnTZXQRhsQ2K6VKAXQXkfoAvhKRzmUsbvltFpGxAHKVUhtEZIg/LzGZZqlt1gxQSh0WkQQAP4nIzjKWrZJttloNPQdAC8PzJACHg1SWqnJMa3ZB+5urTfe27TnaY/fp1ZKIRMIe5v9TSn2pTQ7pbXZQShUAWAZgJEJ7mwcA+I2IZMPeLTpURD5EaG8zlFKHtb+5AL6CvYv4km6z1QJ9PYBUEWktIlEAxgP4JshlCrRvANypPb4TwDzD9PEiEi0irQGkAlinNePOiEhf7Wz4HYbXVCta+d4FsEMp9aphVihvc7xWM4eI1AJwFYCdCOFtVko9qZRKUkolw/4dXaKUuh0hvM0iUkdEYh2PAYwAkIlLvc3BPjNcgTPJo2EfHbEPwFPBLk8lt+VjAEcAFMN+ZJ4IoBGAxQD2aH8bGpZ/StvuXTCc+QaQpu08+wC8Du0K4Or2D8BA2JuPWwBkaP9Gh/g2dwWwSdvmTABTtekhu81u2z8EzlEuIbvNsI+826z92+bIpku9zbz0n4goRFity4WIiLxgoBMRhQgGOhFRiGCgExGFCAY6EVGIYKATEYUIBjoRUYj4f3fBaaT8x+ppAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQsklEQVR4nO3df4hlZ33H8ffH/CAaDdltZpetcR2lSzQIJumQxgak7TZtNOKGYoqCdiuRpVBF24Ks/iP9b6FFtCCWJdGOGNOmMSFLtNYwKlKw0c0PNXGTrsYYU9fdMdbGKlSi3/4xJ3Wcndm5c3/Nfea+XzCce865P77PvXc/+9znPufcVBWSpPY8Z7MLkCT1xwCXpEYZ4JLUKANckhplgEtSo84e54NddNFFNTs7O86HlKTm3XfffT+oqpmV28ca4LOzsxw9enScDylJzUvyndW2O4QiSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsA1NWYPfmqzS5CGygCXpEYZ4JLUKANckhplgEtSo9YN8CSXJHlw2d/TSd6VZHuSe5Ic75bbxlGwJGnJugFeVY9W1WVVdRnwm8BPgTuBg8BCVe0BFrp1SdKYbHQIZS/wrar6DrAPmO+2zwPXD7EuSdI6NhrgbwRu7S7vrKoTAN1yx2o3SHIgydEkRxcXF/uvVJL0K3oO8CTnAq8H/nkjD1BVh6tqrqrmZmZO+0k3SVKfNtIDfw1wf1Wd7NZPJtkF0C1PDbs4SdLaNhLgb+KXwycAR4D93eX9wF3DKkqStL6eAjzJ84BrgDuWbT4EXJPkeLfv0PDLkySt5exerlRVPwV+bcW2p1ialSJJ2gQeiSlJjTLAJalRBrgkNcoAV9NmD37KH2rQ1DLAJalRBrgkNcoA16Zw2EManAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuLRBngFRk8IAl6RG9fqjxhcmuT3JI0mOJXlVku1J7klyvFtuG3WxkqRf6rUH/kHgM1X1MuCVwDHgILBQVXuAhW5dkjQm6wZ4kguAVwM3A1TVz6rqR8A+YL672jxw/WhKlCStppce+EuBReCjSR5IclOS84GdVXUCoFvuGGGdkqQVegnws4ErgA9X1eXAT9jAcEmSA0mOJjm6uLjYZ5mSpJV6CfAngSer6t5u/XaWAv1kkl0A3fLUajeuqsNVNVdVczMzM8OoWZLEUu/6jKrq+0m+m+SSqnoU2At8o/vbDxzqlneNtFJNheXzqx8/dN0mViJNvnUDvPMO4JYk5wKPAW9lqfd+W5IbgSeAG0ZToiRpNT0FeFU9CMytsmvvUKuRJPXMIzElqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrg0Yv6GpkbFAJekRhngktQoA1wjMY3DBtPYZm0uA1ySGtXr+cClqWfvWpPGHrgkNcoAl6RGGeCS1CgDXJIaZYBLUqN6moWS5HHgx8DPgWeqai7JduCfgFngceCPq+q/RlOmJGmljfTAf7eqLquqZ3+d/iCwUFV7gIVuXZI0JoMMoewD5rvL88D1A1cjSepZrwFewGeT3JfkQLdtZ1WdAOiWO1a7YZIDSY4mObq4uDh4xZIkoPcjMa+uqu8l2QHck+SRXh+gqg4DhwHm5uaqjxolSavoqQdeVd/rlqeAO4ErgZNJdgF0y1OjKlKSdLp1AzzJ+Ule8Oxl4A+Ah4AjwP7uavuBu0ZVpCTpdL0MoewE7kzy7PU/UVWfSfIV4LYkNwJPADeMrkxJ0krrBnhVPQa8cpXtTwF7R1GUJGl9Hokp4Y8xqE0GuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLApWU2MhfceePabAa4JDXKAJekRhng0grDPKzeYRaNkgEuSY0ywCWpUQa4JDXKAJekRhngktSoXn+VXhqpaZmt8Ww7Hz903SZXoq3AHrgkNarnAE9yVpIHktzdrW9Pck+S491y2+jK1FbR7xzrjdxuo/c/Lb1/bT0b6YG/Ezi2bP0gsFBVe4CFbl2SNCY9BXiSi4HrgJuWbd4HzHeX54Hrh1qZJOmMeu2BfwB4N/CLZdt2VtUJgG65Y7UbJjmQ5GiSo4uLi4PUKm1Zwzx8X9Nj3QBP8jrgVFXd188DVNXhqpqrqrmZmZl+7kKStIpephFeDbw+yWuB84ALknwcOJlkV1WdSLILODXKQiVJv2rdHnhVvaeqLq6qWeCNwOeq6s3AEWB/d7X9wF0jq1JaxSQNOUxSLZoeg8wDPwRck+Q4cE23Lkkakw0diVlVXwC+0F1+Ctg7/JIkSb3wSExtSc7q0DQwwCWpUQa4JDXKAJca5lDRdDPAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBrQ3qZ8TDKWRHOuFji7BOBAS5JzfJX6aUzWO9X5O0FazPZA5ekRhngktQoh1A00SZliGJS6pCWswcuSY0ywCWpUQ6h6DTrzbyYFoMMmzjkonGwBy5JjTLAJalR6wZ4kvOSfDnJV5M8nOSvu+3bk9yT5Hi33Db6ciUPI98In6etrZce+P8Cv1dVrwQuA65NchVwEFioqj3AQrcuSRqTdQO8lvxPt3pO91fAPmC+2z4PXD+KAiVJq+tpFkqSs4D7gN8APlRV9ybZWVUnAKrqRJIda9z2AHAAYPfu3cOpWhNprY/rzmqRRqOnLzGr6udVdRlwMXBlklf0+gBVdbiq5qpqbmZmps8yJUkrbWgWSlX9CPgCcC1wMskugG55atjFSZLW1ssslJkkF3aXnwv8PvAIcATY311tP3DXiGrUBBtkRoizSU53pufD50sr9TIGvguY78bBnwPcVlV3J/kScFuSG4EngBtGWKckaYV1A7yqvgZcvsr2p4C9oyhKW4c9Rml0PBJTkhplgEtSowxwbbrNGGZxaEdbgQEuSY0ywCWpUf6gg37F8qGFrXIIvMMl2qrsgUtSowxwSWqUAS5tQRs57N4hpnYZ4JLUKANckhrlLJQtaqvMIBmlloYOeqnV13z62AOXpEYZ4JLUKANcY9PSkMUkG8bz6GuxNRjgktQoA1ySGmWAS1KjDHBJalQvv0r/oiSfT3IsycNJ3tlt357kniTHu+W20ZerYZqGXzlv8ccipuF10XD00gN/Bvirqno5cBXw50kuBQ4CC1W1B1jo1iVJY7JugFfViaq6v7v8Y+AY8EJgHzDfXW0euH5ENUqSVrGhQ+mTzAKXA/cCO6vqBCyFfJIda9zmAHAAYPfu3QMVK53JRg43l7aCnr/ETPJ84JPAu6rq6V5vV1WHq2ququZmZmb6qVGStIqeAjzJOSyF9y1VdUe3+WSSXd3+XcCp0ZQoSVrNukMoSQLcDByrqvcv23UE2A8c6pZ3jaRCSUPjENLW0ssY+NXAW4CvJ3mw2/ZeloL7tiQ3Ak8AN4ykQknSqtYN8Kr6NyBr7N473HIkSb3ySEytq4WP3S3U2I9xHdTT7+Ns1ee9FQa4JDXKAJekRhng2jDP1TEdfI0nnwEuSY3yV+mlLWCYveVh97yfvb/HD1031PuVPXBJapYBLkmNMsClTbbWkMVW/BLRL8CHywCXpEYZ4JLUKANcffOj8GRa+bpsxrCFQyXjYYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKc6FI2jBnmEwGe+CS1Kh1AzzJR5KcSvLQsm3bk9yT5Hi33DbaMiVJK/XSA/8H4NoV2w4CC1W1B1jo1tWnYX0c9eAJbQbfd5tn3QCvqi8CP1yxeR8w312eB64fblmSpPX0+yXmzqo6AVBVJ5LsWOuKSQ4ABwB2797d58NpkJPiD+OE+vawpMkz8i8xq+pwVc1V1dzMzMyoH06Spka/AX4yyS6AbnlqeCVJknrRb4AfAfZ3l/cDdw2nHG3EJA1rTFItWt1mvEa+L0arl2mEtwJfAi5J8mSSG4FDwDVJjgPXdOuSpDFa90vMqnrTGrv2DrkWSdIGeCTmhFltTu20fgyd1nZPol5ei0HmgzuXvD8GuCQ1ygCXpEYZ4I0Z5GPmmW7rR1iNm++5wRngktQoA1ySGmWAT4hxf5z0o+vmmtTnv4W6Rj0jpiUGuCQ1ygCXpEb5m5h9GMbpWYfx+JN6f9JG+P7rnz1wSWqUAT5ia/Uu/Bk1tfraDXK4/KD30aJRvs4GuCQ1ygCXpEYZ4Cus93FntTMFTsLHwVEP1Whwrb8Wvb7XR3W6B53OAJekRhngktQoA5zJGQaRxmkrvO/7ndkyqUOhG2WAS1KjDHBJatRAh9InuRb4IHAWcFNVjezX6Xs5fH2j1+n1I9PswU8N7bD5Xj/yDfMxpa1mGMMdwzolxmqZMq5/u333wJOcBXwIeA1wKfCmJJcOqzBJ0pkNMoRyJfDNqnqsqn4G/COwbzhlSZLWk6rq74bJG4Brq+pt3fpbgN+qqrevuN4B4EC3egnwaP/ljt1FwA82u4hNYtunk22fTC+uqpmVGwcZA88q207736CqDgOHB3icTZPkaFXNbXYdm8G22/Zp02LbBxlCeRJ40bL1i4HvDVaOJKlXgwT4V4A9SV6S5FzgjcCR4ZQlSVpP30MoVfVMkrcD/8rSNMKPVNXDQ6tsMjQ59DMktn062faG9P0lpiRpc3kkpiQ1ygCXpEZNZYAnuTbJo0m+meTgKvu3JbkzydeSfDnJK5btuzDJ7UkeSXIsyavGW/1gBmz7XyR5OMlDSW5Nct54qx9Mko8kOZXkoTX2J8nfdc/N15JcsWzfGZ+3Sddv25O8KMnnu/f6w0neOd7KBzfI697tPyvJA0nuHk/FG1BVU/XH0heu3wJeCpwLfBW4dMV1/gZ4X3f5ZcDCsn3zwNu6y+cCF252m8bRduCFwLeB53brtwF/utlt2mD7Xw1cATy0xv7XAv/C0jEOVwH39vq8TfrfAG3fBVzRXX4B8B/T0vZl+/8S+ARw92a3ZeXfNPbAezkFwKXAAkBVPQLMJtmZ5AKW3gw3d/t+VlU/Glvlg+u77d2+s4HnJjkbeB6Nzfuvqi8CPzzDVfYBH6sl/w5cmGQXW+C0Ef22vapOVNX93X38GDjG0n/mzRjgdSfJxcB1wE2jr3TjpjHAXwh8d9n6k5z+hvwq8EcASa4EXszSgUovBRaBj3YfqW5Kcv7oSx6avtteVf8J/C3wBHAC+O+q+uzIKx6vtZ6fXp631q3bxiSzwOXAveMrayzO1PYPAO8GfjHmmnoyjQHeyykADgHbkjwIvAN4AHiGpR7oFcCHq+py4CdAS+Ohfbc9yTaWeiovAX4dOD/Jm0dY62ZY6/np6bQRjTtjG5M8H/gk8K6qenpsVY3Hqm1P8jrgVFXdN+6CejXQ+cAbte4pALo36Fth6QsOlsZ+v83SsMGTVfVsD+R22grwQdr+h8C3q2qx23cH8NvAx0df9tis9fycu8b2rWTN90aSc1gK71uq6o5NqG3U1mr7G4DXJ3ktcB5wQZKPV9XEdFymsQe+7ikAupkm53arbwO+WFVPV9X3ge8muaTbtxf4xrgKH4K+287S0MlVSZ7XBftelsZDt5IjwJ90sxKuYmmY6ATTcdqIVdvevdY3A8eq6v2bW+LIrNr2qnpPVV1cVbMsveafm6TwhinsgdcapwBI8mfd/r8HXg58LMnPWQroG5fdxTuAW7p/yI/R9VZbMEjbq+reJLcD97M0nPQAjR16nORW4HeAi5I8CbwPOAf+v+2fZmlGwjeBn9K9tms9b2NvwAD6bTtwNfAW4OvdsBrAe6vq02MrfkADtH3ieSi9JDVqGodQJGlLMMAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSo/4PZiHUKjh3dY8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 9.19289491 10.45883096 10.80876574]]\n"
     ]
    }
   ],
   "source": [
    "plt.hist(w_b[2000:,4],bins=200)\n",
    "plt.show()\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9986154788807169"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w_b[:,4].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPo0lEQVR4nO3df6xfdX3H8eeLguBPhFB+WMrKXMcEM2G7YW4kiwtuMCcWk7mUbIRsuroENlzMNnB/6P5oYqLiTDZNKjBJRBlBVOqciszEaAQsPyKU2lmBwYVCO5mAv4C27/1xT83X9t7eH9/zvd/eT5+P5OZ+v59zzve8T9r7up/7Oed8TqoKSVJbDht3AZKk/hnuktQgw12SGmS4S1KDDHdJatDh4y4A4LjjjqtVq1aNuwxJWlLuuuuu/62q5dMtOyjCfdWqVWzatGncZUjSkpLkf2Za5rCMJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ16KC4Q1U6JFxwwfTtGzcubh06JNhzl6QGGe6S1CDDXZIaNGu4J1mZ5GtJtiTZnOTyrv19SR5Lcm/39aaBba5Msi3J1iTnjfIAJEn7m8sJ1V3Au6vq7iQvB+5Kcmu37MNV9cHBlZOcDqwFzgBeBXw1ya9W1e4+C5ckzWzWnntVba+qu7vXzwJbgBUH2GQNcENVPVdVDwHbgLP7KFaSNDfzGnNPsgo4C7ija7osyXeSXJvkmK5tBfDowGaTTPPLIMm6JJuSbNq5c+f8K5ckzWjO4Z7kZcBngHdV1TPAx4BXA2cC24EP7V11ms1rv4aqDVU1UVUTy5dP+5QoSdICzekmpiRHMBXs11fVzQBV9eTA8o8DX+jeTgIrBzY/GXi8l2qlFnlzk0ZgLlfLBLgG2FJVVw20nzSw2luB+7vXtwBrkxyZ5FRgNXBnfyVLkmYzl577OcDFwH1J7u3a3gNclORMpoZcHgbeCVBVm5PcCDzA1JU2l3qljCQtrlnDvaq+wfTj6F88wDbrgfVD1CVJGoJ3qEpSgwx3SWqQU/5KfZvp6hdpEdlzl6QGGe6S1CDDXZIaZLhLUoM8oSrNxukBtATZc5ekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhrk3DLSXj5kQw2x5y5JDTLcJalBDstIC+Uwjg5i9twlqUGGuyQ1yHCXpAYZ7pLUIE+oSkuNz3TVHNhzl6QGGe6S1CCHZXTo8fp0HQJm7bknWZnka0m2JNmc5PKu/dgktyb5Xvf9mIFtrkyyLcnWJOeN8gAkSfuby7DMLuDdVfUa4PXApUlOB64Abquq1cBt3Xu6ZWuBM4DzgY8mWTaK4iVJ05s13Ktqe1Xd3b1+FtgCrADWANd1q10HXNi9XgPcUFXPVdVDwDbg7J7rliQdwLzG3JOsAs4C7gBOqKrtMPULIMnx3WorgNsHNpvs2vb9rHXAOoBTTjll3oVL2oeXSGrAnK+WSfIy4DPAu6rqmQOtOk1b7ddQtaGqJqpqYvny5XMtQ5I0B3MK9yRHMBXs11fVzV3zk0lO6pafBOzo2ieBlQObnww83k+5kqS5mHVYJkmAa4AtVXXVwKJbgEuA93ffPz/Q/qkkVwGvAlYDd/ZZtA5BDjlI8zKXMfdzgIuB+5Lc27W9h6lQvzHJ24FHgLcBVNXmJDcCDzB1pc2lVbW778IlSTObNdyr6htMP44OcO4M26wH1g9RlyRpCN6hKh2svJNWQzDc1SaDUYc4Jw6TpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CCnH9DS5jQD0rTsuUtSgwx3SWqQ4S5JDTLcJalBnlCVDmU+m7ZZ9twlqUH23HVw8dJGqRf23CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAbNGu5Jrk2yI8n9A23vS/JYknu7rzcNLLsyybYkW5OcN6rCJUkzm0vP/RPA+dO0f7iqzuy+vgiQ5HRgLXBGt81Hkyzrq1hJ0tzMGu5V9XXgqTl+3hrghqp6rqoeArYBZw9RnyRpAYYZc78syXe6YZtjurYVwKMD60x2bftJsi7JpiSbdu7cOUQZkqR9LTTcPwa8GjgT2A58qGvPNOvWdB9QVRuqaqKqJpYvX77AMiRJ01nQwzqq6sm9r5N8HPhC93YSWDmw6snA4wuuTtLwfADKIWlBPfckJw28fSuw90qaW4C1SY5MciqwGrhzuBIlSfM1a889yaeBNwDHJZkE3gu8IcmZTA25PAy8E6CqNie5EXgA2AVcWlW7R1K5JGlGs4Z7VV00TfM1B1h/PbB+mKIkScPxAdmS9jfTOP3GjYtbhxbM6QckqUH23DUeXsEhjZQ9d0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yInDNFpOECaNhT13SWqQ4S5JDTLcJalBhrskNchwl6QGebWMpOH5QO2DjuEuae68tHXJcFhGkhpkz13S6DhcMzb23CWpQYa7JDXIcJekBhnuktQgT6iqH14iJx1U7LlLUoNmDfck1ybZkeT+gbZjk9ya5Hvd92MGll2ZZFuSrUnOG1XhkqSZzaXn/gng/H3argBuq6rVwG3de5KcDqwFzui2+WiSZb1VK0mak1nDvaq+Djy1T/Ma4Lru9XXAhQPtN1TVc1X1ELANOLufUiVJc7XQE6onVNV2gKranuT4rn0FcPvAepNd236SrAPWAZxyyikLLEOLzhOn0pLQ9wnVTNNW061YVRuqaqKqJpYvX95zGZJ0aFtouD+Z5CSA7vuOrn0SWDmw3snA4wsvT5K0EAsN91uAS7rXlwCfH2hfm+TIJKcCq4E7hytRkjRfs465J/k08AbguCSTwHuB9wM3Jnk78AjwNoCq2pzkRuABYBdwaVXtHlHtkqQZzBruVXXRDIvOnWH99cD6YYqSJA3HO1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg3xYh6TFN9McRRs3Lm4dDbPnLkkNMtwlqUGGuyQ1yDF3Tc9526UlzZ67JDXIcJekBhnuktQgw12SGuQJVUkHD29u6o09d0lqkOEuSQ1yWOZQ5/XsUpPsuUtSg+y5HyrsoUuHFHvuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho01PQDSR4GngV2A7uqaiLJscC/A6uAh4E/qar/G65MSdJ89NFz/72qOrOqJrr3VwC3VdVq4LbuvSRpEY1iWGYNcF33+jrgwhHsQ5J0AMOGewFfSXJXknVd2wlVtR2g+378dBsmWZdkU5JNO3fuHLIMSdKgYaf8PaeqHk9yPHBrku/OdcOq2gBsAJiYmKgh65AkDRgq3Kvq8e77jiSfBc4GnkxyUlVtT3ISsKOHOjVXztsuiSGGZZK8NMnL974G/gC4H7gFuKRb7RLg88MWKUman2F67icAn02y93M+VVVfSvJt4MYkbwceAd42fJmSpPlYcLhX1YPA66Zp/wFw7jBFSZKG4x2qktQgH5C9FHnSVNIsDHdJS9d8OzobN46mjoOQwzKS1CB77pIOfg5Fzps9d0lqkOEuSQ1yWOZg5p+ikhbInrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBjn9wMHAaQakxbGQn7UlOge84b6YDHFJi8RhGUlqkD33UbCHLmnM7LlLUoMMd0lqkMMyktSnmYZlF/mqG3vuktQge+7D8MSp1L6+euKL3KO35y5JDTLcJalBDsvMhcMvkvZ1kOeCPXdJatDIeu5Jzgc+AiwDrq6q9/e9j8/d8xgf+PJWHvvhT1mWsLuKFa98MX933mlceNYKuOACnnj6Z2x98lle2L1n2s84YtlhnHbCyznx6KN+3vatB3/Aj5/b9fP3Rx5+GEn42Qu7D1jPYQCBPdXH0UlaapYdFo5+8RE89ePn57zNO674DwDOefWxXP+Xv91bLSMJ9yTLgH8Ffh+YBL6d5JaqeqCvfXzunse48ub7+OkLu7n6pn/6hWXLrg1PnPgKAB7Y/gx7aua0fWH3Hh7Y/gwAJx591H7BDvDcrul/MexrD4DBLh2ydu+peQX7oG9+/yn+9OPf6i3gRzUsczawraoerKrngRuANX3u4ANf3spPZ+hJ795TbNv5I7bt/NEBg32vPTW1PrBfsEvSYvnm95/q7bNScwi/eX9o8sfA+VX1ju79xcBvVdVlA+usA9Z1b08Dts5nHy868Vd+s6dyf+75J7bdtZDP3f2Tp1n2kqP7LmfsPK6lp9VjO5SO6/kntt01j4/4papaPt2CUY25Z5q2X/gtUlUbgA0j2v+iSrJp19M7JsZdR988rqWn1WPzuOZvVMMyk8DKgfcnA4+PaF+SpH2MKty/DaxOcmqSFwFrgVtGtC9J0j5GMixTVbuSXAZ8malLIa+tqs2j2NdBoonhpWl4XEtPq8fmcc3TSE6oSpLGyztUJalBhrskNchw70mSDyT5bpLvJPlskleOu6ZhJDk/ydYk25JcMe56+pBkZZKvJdmSZHOSy8ddU5+SLEtyT5IvjLuWPiV5ZZKbup+vLUn6u0d/jJL8bff/8P4kn05y1OxbzZ3h3p9bgddW1a8D/w1cOeZ6Fmxg+og/BE4HLkpy+nir6sUu4N1V9Rrg9cCljRzXXpcDW8ZdxAh8BPhSVf0a8DoaOMYkK4C/ASaq6rVMXXiyts99GO49qaqvVNXeuQtuZ+ra/qVq5NNHjENVba+qu7vXzzIVEivGW1U/kpwM/BFw9bhr6VOSVwC/C1wDUFXPV9UPx1pUfw4HXpzkcOAl9HwvkOE+Gn8B/Oe4ixjCCuDRgfeTNBKCeyVZBZwF3DHmUvryz8Df081f15BfBnYC/9YNOV2d5KXjLmpYVfUY8EHgEWA78HRVfaXPfRju85Dkq9342L5fawbW+Uem/vy/fnyVDm3W6SOWsiQvAz4DvKuqnhl3PcNK8mZgR1XNZ06SpeJw4DeAj1XVWcCPgSV/DijJMUz9NXwq8CrgpUn+rM99+CSmeaiqNx5oeZJLgDcD59bSvoGg2ekjkhzBVLBfX1U3j7uenpwDvCXJm4CjgFck+WRV9RoWYzIJTFbV3r+wbqKBcAfeCDxUVTsBktwM/A7wyb52YM+9J93DSf4BeEtV/WTc9QypyekjkoSpsdstVXXVuOvpS1VdWVUnV9Uqpv6t/quRYKeqngAeTXJa13Qu0NtzIcboEeD1SV7S/b88l55PFNtz78+/AEcCt079W3F7Vf3VeEtamIanjzgHuBi4L8m9Xdt7quqL4ytJc/DXwPVdR+NB4M/HXM/QquqOJDcBdzM1jHsPPU9F4PQDktQgh2UkqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQ/wMtwwFWx8ikmQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import random\n",
    "import math\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "def norm_dist_prob(theta):\n",
    "    y = norm.pdf(theta, loc=3, scale=2)\n",
    "    return y\n",
    "\n",
    "T = 5000\n",
    "pi = [0 for i in range(T)]\n",
    "sigma = 1\n",
    "t = 0\n",
    "while t < T-1:\n",
    "    t = t + 1\n",
    "    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)\n",
    "    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))\n",
    "\n",
    "    u = random.uniform(0, 1)\n",
    "    if u < alpha:\n",
    "        pi[t] = pi_star[0]\n",
    "    else:\n",
    "        pi[t] = pi[t - 1]\n",
    "\n",
    "\n",
    "plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))\n",
    "num_bins = 50\n",
    "plt.hist(pi, num_bins, facecolor='red', alpha=0.7)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "d83df975d3fd7e9cd2b9393f36d258bee244a85bb4d399c6c89ec18ff8dd5f85"
  },
  "kernelspec": {
   "display_name": "Python 3.8.3 ('data')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}